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A method of numerically estimating dynamic Green’s functions using the finite element 
method is proposed. These Green’s functions are accurate in a limited frequency range 
dependent on the mesh size used to generate them. This range can often match or exceed 
the frequency sensitivity of the traditional acoustic emission sensors. An algorithm is also 
developed to characterize an acoustic emission source by obtaining information about its 
strength and temporal dependence. This information can then be used to reproduce the 
source in a finite element model for further analysis. Numerical examples are presented 
that demonstrate the ability of the band-limited Green’s functions approach to determine 
the moment tensor coefficients of several reference signals to within seven percent, as well 
as accurately reproduce the source-time function. 


Nomenclature 


G Dynamic Green’s function 

x,f Location vectors 

t Time 

Cijki Elastic stiffness tensor 
u.i Displacement vector 
M Seismic moment tensor 
?n,; Moment tensor coefficient 
* Convolution operator 
A Bulk modulus 

p Shear modulus 


I. Introduction 

Acoustic emissions are transient elastic waves resulting from the sudden release of strain energy at the 
onset and/or progression of material damage or degradation. These elastic waves propagate through the 
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medium where they can be detected by sensors mounted remotely on the surface. Acoustic emission (AE) 
monitoring can be used as a viable nondestructive method for the real time monitoring of damage progression 
in a structure while in service. AE differs from other nondestructive evaluation (NDE) methods in that it 
monitors generally irreversible processes and therefore can only be used to detect the evolution of damage 
as it occurs. The advantage of AE is that it only requires a small set of passive transducers to monitor a 
structure. 1 

AE analysis has typically been considered a qualitative exercise, where correlations between experimental 
results and AE signal characteristics are observed and used in pattern recognition schemes to draw conclusions 
during future AE monitoring applications. 1 The limitation with this approach is that there is no underlying 
physics identified between the experimental results and the AE signal characteristics. Different sources 
producing AE signals with similar characteristics may lead to unsubstantiated conclusions being drawn. To 
overcome this shortcoming, methods of quantitative AE evaluation have been developed to study the details 
of the AE source and resulting waves. These methods will allow more quantitative conclusions to be drawn 
about the correlation between damage and the resulting AE signals as well as providing for an improved 
ability to distinguish between sources that share similar characteristics. 2,3 

II. Theory 

II. A. Theoretical Representation of an Acoustic Emission Source 

To extract quantitative data from an acoustic emission event, the source of the emission must be able to 
be mathematically modeled in some fashion. Scruby has shown that an acoustic emission source can be 
represented as a point source with two main assumptions. 2 First, it is assumed the source is internal and 
self-equilibrating. This means that every force or couple that makes up the source must be a force dipole or 
double couple. Second, all forces that make up the source occur simultaneously, such that they all share the 
same time history. The displacement response on the surface of the body due to the point source can then 
be calculated as follows. 

Consider a force, applied in the jth direction to a body at position £ as shown in figure 1. The 

displacement at location x in the ith direction, iq (x, t ) , can be calculated as 


Ui(x,t) = (1) 

where * denotes convolution, vectors are represented with an arrow, and Gij(x,t, £) is the Green’s function, 
defined as the displacement response at x in the ith direction due to an impulse load at £ in the jth direction. 3 



Figure 1. Point source and sensor locations in a general elastic body. 


In order to meet the assumptions of an internal, self-equilibrating source, the force, hj(£,f) is replaced 
with a moment tensor, M, containing the strengths of the system of dipole forces representing the AE source, 
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the components of which are denoted as nijk{£,)- The equivalent force dipole or couple for each moment tensor 
component is shown in figure 2. All of the forces that make up M act with the same temporal behavior, 
/(<). Thus, Eq. (1) becomes 


3 3 

*(*.*) = EE ) * G ij>k (x, t, |) (2) 

j—i fc= l 

where G ih k{x. t, <f) is the Green’s function of a dipole source, defined as the impulse displacement response 
in the ?’th direction at x due to a dipole impulse force in the jth direction at location £, where the impulse 
forces that make up the dipole are separated in the fcth direction. 



Figure 2. Force dipoles and double couples corresponding to the moment tensor. 


The values of the moment tensor coefficients for a displacement discontinuity can be determined for 
homogeneous materials according to the relation 1,3 

Trijj CijkibkdAi (3) 

where Cijki is the elastic stiffness tensor. For an isotropic solid, Cijki is given as 

Gijki = ^bijbfd T /A A/.: A/ -f- Sj t 6 jk ) (4) 

where 5ij is the Kronecker delta = 1 if and only if * = j), and A and /r are Lame’s constants. 

As an example, consider a crack with faces each having surface area dA in the X\X 2 plane that opens in 
the X 3 direction. The equivalent body forces are found by substituting l = k = 3 into Eq. (3) and Eq. (4) 
and the components of the resulting moment tensor are: 


Tfljk 


M 


X5ij + fi(S l 3 S i3 + Sj 3 Sj 3 )bdA 


AO 0 

0 A 0 

0 0 A + 2/z 


bdA 


( 5 ) 

( 6 ) 
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II. B. Finite Element Representation of an Acoustic Emission Source 

The finite element method (FEM) can be used to model acoustic emission wave forms. FEM has several 
distinct advantages over traditional analytical approaches to modeling AE waveforms. Analytical approaches 
are typically limited to very simple geometry, i.e. thin, infinite plates or half spaces. In contrast, FEM can 
handle complex geometries, including boundary reflections. The main disadvantage of using FEM is the 
computational time required. The frequency range that can be modeled is determined by the element size 
used (smaller elements are needed to capture higher frequencies), so there is a trade off that must be made 
between the frequency range of the model and the computational time required. 4 

AE sources cannot be modeled as a point source when using finite elements. Therefore, sources are 
approximated as point sources by body forces that are applied over some finite area. Hamstad, et. al. have 
developed a set of criteria for the convergence of a finite element model of an acoustic emission source based 
on the maximum frequency of interest, f ma x, which is inversely proportional to the minimum wavelength of 
interest, A min . 4 


A 


min 


Ct 


( 7 ) 


where ct is the transverse, or shear, wave velocity in the material. Several relations were developed between 
the minimum wavelength, A mini source size, s, finite element mesh size, MS, and the distance between 
the source and the transducer, D. Convergence was measured by the normalized amplitude of results with 
different parameters. The following criteria were found for convergence of peak amplitudes in modeled signals 
to within 5%, 4 


An 


MS 

D 


> 

2 

(8) 

> 

15 

(9) 

> 

7 

(10) 


In this study, all FEM calculations are made using the commercially available software ABAQUS™. 5 
A comparison between an analytical solution, found in the literature, for a buried dipole in a steel plate 
and the results of a finite element simulation is shown in figure 3. 4 It can be seen that the FEM solution 
accurately captures all of the features of the analytical solution. Therefore, FEM will be used to calculate 
all reference solutions for this study. For the comparison, a vertical force dipole of strength 0.3 N ■ mm. was 
located at the center of an axi-symmetric model of a 25 mm thick steel disk. The sensor was located 50 mm 
from the epicenter as shown in figure 4. FEM parameters for all simulations done in this study are shown 
in table 1. 


Table 1. Model parameters for FEM analysis. 



Mesh (/im.) 

Elements 

Dipole Spacing (pm)* 

E (GPa) 

V 

P (kg/m 3 ) 

Figure 4: 

50 

4 node axi-symmetric 

50 

209 

0.29 

7800 

Figure 5: 

50 

4 node axi-symmetric 

50 

73 

0.30 

2770 

Figure 12: 

50 

8 node brick 

50 

73 

0.30 

2770 

Figure 14: 

100** 

10 node tet 

50 

73 

0.30 

2770 


*A11 sources are 150 pm in size. 50 /. im for each dipole with 50 jirn separating them. 

**Meslr size at source equal to 50 /im, constrained to a maximum size of 100 /im everywhere else. 


II. C. Calculation of the Dynamic Green’s Function 

Central to the characterization of AE events, is the calculation of the dynamic Green’s function between the 
source location and each sensor. Analytical solutions to Green’s functions are difficult if not impossible to 
determine for all but the simplest of geometries, including infinite plates and half-spaces. 2 It is proposed 
that numerical approximations of the dynamic Green’s functions for a source/sensor pair can be calculated 
using FEM. The numerical Green’s functions calculated in this manner are not true Green’s functions. A 
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Figure 3. Comparison between a finite element generated AE signal and an analytically derived signal in a 25 mm 
thick steel disk from Hamstad (1999). 



Figure 4. Schematic of geometry and boundary conditions for axi-symmetric finite element model used for comparison 
with analytical solution. 


true Green’s function will represent the impulse response of the system at all frequencies. The numerical 
approximations of the Green’s function generated by FEM however, are band-limited based on the mesh size 
used to generate them. Therefore, any AE signal modeled using a band-limited Green’s function will not 
converge at higher frequencies. However, the typical AE sensor loses sensitivity in the 1-2 MHz range. 5 6 For 
the quantitative analysis of experimental AE data, the frequency range of the band-limited Green’s function 
only needs to match that of the sensor being used. 

For convenience, a modification is made to Eq. (2). First, it is assumed that the sensors used in all 
experiments are only sensitive to displacement normal to the surface of the body they are monitoring. 
Therefore, the index i in Eq. (2) can be dropped as it is assumed that the displacement, u , and Green’s 
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function, G, always reference the displacement normal to the surface. Second, due to the symmetry of linear 
elasticity, an index transformation can be applied to nijk and Gys, such that, 


mu 

= m 1 

m 2 2 

= m 2 

m .33 

= m 3 

mi2&m 2 i 

= 777.4 


= 7775 

m 2 3&m 3 2 

= m 6 

Gi,i 

= Gi 

G 2 ,2 

= g 2 

G3 , 3 

= g 3 

Gi j 2&G2,1 

= g 4 

g 1j3 &g 3 , 1 

= g 5 

G2, 3 &G 3j 2 

= g 6 


thus, Eq. (2) becomes 

6 

u r = mjf(t)*G r j r=l...n (11) 

l=i 

where n is the number of sensors. 

A finite element model of the specimen is then generated and run six times in order to obtain a set of 
displacements related to each component of the moment tensor. Each run uses a different dipole or double 
couple of a known strength and time dependence located at the position of the AE source. For example, 
for the first run set m\ = 1 and 7772-6 = 0. For the second run, m 2 = 1 and mi = 7?7 3 _6 = 0 and so on. 
Displacement output normal to the surface is recorded at each sensor location. The displacement output 
corresponds to u r and the prescribed dipole strength and time dependence correspond to rrij and /(f) in 
Eq. (11). Then, Gj is found through the deconvolution of from u r . 

For the present study, a least squares approach to deconvolution is used. Let j = 1 and mi = 1. The 
convolution equation representing the displacement response due to a dipole in the x\ direction is given by 


u(t) = f{t) * Gi(t) 

and can be rewritten in matrix form as 3 


Ui 

U 2 


fi 0 0 . 

h h 0 . 

0 

0 


' G u ' 

Gi 2 

Um 


fm fm— 1 fm — 2 

• • 1 


. G U . 


(12) 


(13) 


where m and n are the number of the data points in the displacement and Green’s function vectors. The 
unknown Green’s function Gift) is then found by solving Eq. (13) using the conjugate gradient method.' 
Deconvolution is extremely sensitive to error, especially in the higher frequency range. Therefore, it is 
important that the tolerances are not set too small during the deconvolution process. Some trial and error 
may be required to obtain optimal performance. 8 


II. D. Calculation of the Moment Tensor and Source-Time Dependence 

Once the band-limited Green’s functions have been found for each component of the moment tensor, the 
moment tensor and source time dependence can be recovered from Eq. (11). This is accomplished by making 
an initial guess for f(t) and calculating the values of m,j using the conjugate gradient method as follows. 
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Let Wj represent the convolution of each Green’s function with the initial guess of the source time function 


Wj = f(t)* Gj 


(14) 


Convolution in the time domain is equivalent to multiplication in the frequency domain, so to calculate Wj, 
f(t ) and Gj are transformed into the frequency domain using the Fourier transform, multiplied together, 
and then returned to the time domain via the inverse Fourier transform. Substituting Eq. (14) into Eq. (11) 
gives 


6 

u r = ^2 m jWj r=l...n (15) 

j'=i 

Let [u(i)] denote the time integral of function u(t) such that 

[u(t)] = f u(t)dt (16) 

Jo 

Integrating both sides of Eq. (15) yields 

[i( r ] = mi[w(] + m 2 [w 2 ] + m3 [11)3] + 7714(^4] + m 5 [wl] + me [wg] r = l...n (17) 

which can be written in matrix form as 


M 1 ] 


R] 

R] • 

■ • R] 


TOi 


U 2 ] 

= 

k] 

Rj • 

■ • R] 


m 2 

(18) 

U n } 


. K] 

K] • 

•• K] . 


. TO 6 



Eq. (18) can then be solved using the conjugate gradient method. The resulting values of mj are then used 
to deconvolve a refined approximation of /(f) from Eq. (11) using the same method applied in Eq. (13) that 
was used to find the band-limited Green’s functions. The process is repeated until the values of m, and /(f) 
converge. 


III. Numerical Examples 


III. A. Estimation of the Dynamic Green’s Function 

A band-limited Green’s function was calculated for a buried dipole at the mid-plane of a 10 mm thick 
aluminum disk of radius 50 mm. The sensor was located 25 mm from the center of the plate as shown in 
figure 5. The Green’s function was then used with Eq. (11) to calculate the response from 3 different dipole 
sources. The moment tensor (or dipole strength) used for all 3 sources was 


M = 


0 0 0 

0 0.3 0 

0 0 0 


(19) 


with the plate laying in the X\X$ plane. The source-time dependence functions for each source are shown 
in figure 6. These functions were chosen arbitrarily to test the Green’s function response to source-time 
functions with differing behavior. In case 1, the dipole forces reach their maximum value after 1.0 /us and 
then persist at a constant value. Case 2 models an impulse with a duration of 1.0 /us. The dipole reaches its 
maximum strength at 0.5 /us and then returns to 0 at 1.0 /us. In case 3, the dipole reaches full strength in 
0.5 /us and then reverses direction and reaches full strength again at 1.5 /us and remains constant. 

To check the accuracy of the numerically estimated Green’s function, the three AE signals calculated 
using Eq. (11) are then compared with reference solutions calculated via direct FEM simulations in figures 7 
(case 1), 8 (case 2) and 9 (case 3). The results show good agreement between the signals calculated with the 
band-limited Green’s function and the signals generated directly by FEM. The maximum difference between 
the two signals is only 0.4%. This demonstrates that the FEM generated Green’s function accurately captures 
the impulse response of the plate. 
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Figure 5. Schematic of geometry and boundary conditions for axi-symmetric finite element model used for numerical 
estimation of the Green’s function for a buried dipole source. 


Case 1 Case 2 Case 3 



Figure 6. Three source-time functions used to simulate buried dipoles in an aluminum plate. 


III.B. Source Characterization 

In order to test the effectiveness of the source characterization algorithm described in section II. D, five AE 
sources were simulated using FEM. The outputs from these models were then analyzed and the moment 
tensor components and source-time functions were determined. First, the three dipoles simulated in figures 
7-8 were analyzed. Next, a general 3-D source was modeled with a prescribed moment tensor and source-time 
function. Finally, a penny-shaped crack was modeled in an aluminum plate. The moment tensor and source 
time function were unknown for this case. 
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Figure 7. Case 1: Comparison between a finite element generated AE signal and the same signal calculated using the 
band-limited Green’s function. 



Figure 8. Case 2: Comparison between a finite element generated AE signal and the same signal calculated using the 
band-limited Green’s function. 


III.B.l. Buried Dipoles 

The three dipoles were analyzed and the resulting calculated dipole strengths are shown in table 2. The 
calculated source-time functions are shown in figure 10. The source-time functions found from the deconvo- 
lution procedure do not all have a peak value of 1. The source-time function found for case 2 has a peak of 
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Figure 9. Case 3: Comparison between a finite element generated AE signal and the same signal calculated using the 
band-limited Green’s function. 


approximately 0.18 and case 3 peaks at 1.6. In both of these cases however, the recovered dipole strengths 
are such that when multiplied by the recovered source-time functions, the resulting peaks are close to the 
prescribed value of 0.3. The calculated source-time functions normalized such that their peak values equal 
1 are shown in figure 11. The corresponding normalized dipole strengths are then shown in the last row 
of table 2. The results match the prescribed values to within 7%. It can also be seen that the recovered 
source-time functions do not settle out at 1 in case 1 and 3, or 0 in case 2, but rather oscillate around these 
values. This is an artifact of the deconvolution procedure. 

Table 2. Calculated dipole strengths 



Case 1 

Case 2 

Case 3 

Prescribed: 

0.30 

0.30 

0.30 

Calculated: 

0.31 

2.00 

0.17 

Normalized: 

0.31 

0.28 

0.28 


% Difference: 3.3 6.6 6.6 


III.B.2. Three Dimensional Source 

A three-dimensional source was modeled in an aluminum plate as shown in figure 12. The source and sensor 
locations are shown in table 3. The prescribed moment tensor coefficients are listed in the first row of table 4. 
The calculated source-time function was badly distorted due the deconvolution process, however, an estimate 
of the rise time (the time it takes the source function to reach an amplitude of 1) could be made. From this 
estimate, a refined initial guess was made that yielded a more accurate result. The prescribed source-time 
function along with the first pass and refined results are shown in figure 13. The source-time function as 
calculated on the first attempt (shown in cyan) reveals the general shape of the true source-time function: 
a rise time of around 0.9 /xs and constant force oscillating around an amplitude of 1 for t > t r i se . Using 
this information, a refined initial guess for the source-time function was made, yielding the results shown 
in red which closely matches the prescribed source-time function. This demonstrates the sensitivity of the 
deconvolution algorithm to the initial guess supplied. The values of the calculated moment tensor coefficients 
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Case 1 Case 2 Case 3 



Figure 10. Comparison between prescribed and recovered source-time functions for three dipole sources. 


Case 1 Case 2 Case 3 



Figure 11. Comparison between prescribed and recovered source-time functions for three dipole sources after normal- 
ization of the recovered functions. 


from the first pass and refined initial guess are shown in table 4 along with the % difference between the 
prescribed values and the final calculated values. The moment tensor coefficients were recovered to within 
3%. 
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Figure 12. Schematic of geometry and boundary conditions for finite element model used for a general 3D source. 


Table 3. Source and sensor locations for 3D source model 



a,’i (mm) 

x 2 (mm) 

x 3 (mm 

Source: 

1.025 

2.525 

0.675 

Sensor 1: 

-15.0 

0.0 

1.0 

Sensor 2: 

0.0 

15.0 

1.0 

Sensor 3: 

15.0 

0.0 

1.0 

Sensor 4: 

0.0 

-15.0 

1.0 

Sensor 5: 

-15.0 

0.0 

0.0 

Sensor 6: 

0.0 

15.0 

0.0 

Sensor 7: 

15 

0.0 

0.0 

Sensor 8: 

0.0 

-15.0 

0.0 

Sensor 9: 

-10.6 

10.6 

1.0 

Sensor 10 

-10.6 

10.6 

0.0 

Sensor 11 

10.6 

10.6 

1.0 

Sensor 12 

10.6 

10.6 

0.0 

Sensor 13 

10.6 

-10.6 

1.0 

Sensor 14 

10.6 

-10.6 

0.0 

Sensor 15 

-10.6 

-10.6 

1.0 

Sensor 16 

-10.6 

-10.6 

0.0 


Table 4. Calculated moment tensor components for a three dimensional AE source 



m 1 

m 2 

m 3 

777,4 

TO 5 

m 6 

Prescribed: 

0.40 

1.35 

0.75 

0.00 

0.00 

0.00 

Calculated: 

0.33 

1.30 

0.58 

0.00 

0.00 

0.00 

Refined: 

0.39 

1.35 

0.73 

0.00 

0.00 

0.00 

% Difference: 

1.65 

0.26 

2.83 

0.00 

0.00 

0.00 


III.B. 3 . Penny Shaped Crack 

A crack was modeled in a small aluminum plate, 2 mm x 2 mm x 1 mm, as shown in figure 14. The faces 
of the crack were held together while the specimen was loaded with a displacement boundary condition of 
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Figure 13. Comparison between an prescribed and recovered source-time function from a three dimensional AE source. 
The impact of an improved initial guess on the accuracy of the solution is demonstrated by the difference between the 
initial calculation shown in cyan and the refined calculation shown in red. 


0.01 mm. After loading, the faces were released and the surface displacements from the sudden opening 
were recorded at 6 locations (listed in table 5). The signals were then analyzed. The recovered source-time 
function is shown in figure 15. The crack was found to have a rise time of approximately 0.3 /is and the 
recovered moment tensor was 


M = 


— 1.73el0 
5.75el0 
0.02el0 


5.75el0 
— 3.62el0 
— 0.20el0 


0.02el0 
— 0.20el0 
1.21el0 


The true source-time function and moment tensor are unknown for this source, therefore, in order to de- 
termine the accuracy of the solution, the recovered moment tensor and source-time function were used to 
reconstruct the signal using Eq. (11). The first 10 /i.s of the reference signal and the reconstruction are shown 
in figure 16. The percent difference between the samples in each signal was calculated. There was an average 
difference of 0.5% with a maximum difference of 6.36%. The performance of the source characterization 
algorithm for the completely unknown source was comparable to the performance for the dipole and 3D 
source. 


Table 5. Source and sensor locations for crack model 

x± (mm) X 2 (mm) X 3 (mm) 


Crack: 

0.0 

0.0 

0.5 

Sensor 1: 

0.0 

0.5 

1.0 

Sensor 2: 

0.0 

-0.5 

1.0 

Sensor 3: 

0.5 

0.0 

0.0 

Sensor 4: 

-0.5 

0.0 

0.0 

Sensor 5: 

1.0 

0.0 

0.5 

Sensor 6: 

-1.0 

0.0 

0.5 
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Figure 14. Schematic of geometry and boundary conditions for finite element model used to model a penny-shaped 
crack. 


IV. Summary 

In this study, a method for estimating band-limited Green’s functions using the finite element method is 
proposed. FEM is shown to be capable of accurately reproducing the displacement response to AE sources. 
Numerical examples were used to demonstrate the accuracy of the band-limited Green’s functions. One FEM 
generated Green’s function was used to calculate the response of 3 dipole sources with different temporal 
dependencies. The calculated responses were compared with FEM generated reference solutions. 

An algorithm was developed to use the band-limited Green’s functions to recover an unknown moment 
tensor and source-time function from an AE signal. Numerical examples were presented in which the source- 
time functions and moment tensors were extracted from FEM generated reference signals. In all cases, the 
rise times of the sources were precisely recovered. At times greater than the rise time, the recovered solutions 
tended to oscillate around the final amplitude of the source due to a build up of high frequency error during 
the deconvolution process. The moment tensor coefficients of the reference signals were recovered within 7 % 
of the prescribed values. A moment tensor and source function were also recovered from a simulated crack. 
Future work will include experimental studies to improve the use of quantitative AE analysis techniques in 
structural health monitoring systems. 
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Normalized Displacement 



Figure 15. Recovered source-time function from the crack model. 



Figure 16. Comparison between the AE signal from the crack model and the reconstructed signal. 
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